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We study noisy heterogeneous diffusion processes with a position dependent diffusivity of the form 
D(x) ~ D(f\x\°‘ in the presence of annealed and quenched disorder of the environment, corresponding 
to an effective variation of the exponent a in time and space. In the case of annealed disorder, 
for which effectively a = aft) we show how the long time scaling of the ensemble mean squared 
displacement (MSD) and the amplitnde variation of individnal realizations of the time averaged MSD 
are affected by the disorder strength. For the case of quenched disorder, the long time behavior 
becomes effectively Brownian after a number of jumps between the domains of a stratified medium. 

In the latter situation the averages are taken over both an ensemble of particles and different 
realizations of the disorder. As physical observables we analyze in detail the ensemble and time 
averaged MSDs, the ergodicity breaking parameter, and higher order moments of the time averages. 

PACS numbers: 05.40.-a 


I. INTRODUCTION 

The motion of individual molecules and submicron 
tracer particles of different sizes in the cytoplasm of living 
biological cells [T] , in artificially crowded environments in 
vitro [ 2 ] , in glass-like systems [ 3 ] , or in large scale in silico 
studies of membrane structures [4] was shown to follow 
the anomalous diffusion law 

( 1 ) 

with the subdiffusive diffusion exponent mostly in the 
range /3 = 0.4...0.9 0 E]. A number of mathemat¬ 
ical models of different kinds were proposed to unveil the 
properties of anomalous diffusion phenomena embodied 
in the mean squared displacement (MSD) in Eq. 0 Cl¬ 
in most of these models the properties of the stochastic 
process are homogeneous in space. Especially for smal¬ 
ler tracers, which may cover longer distances within the 
measurement time, or for techniques allowing for full 
maps of local diffusivities, it turns out that the diffusion 
coefficient becomes a function of the local tracer position. 
For both eukaryotic [5] and prokaryotic [3] cells such local 
diffusivity maps indeed show significant variations. The 
motion of tracer particles through space may also be im¬ 
peded by caging effects when the size of the particle is 
comparable to the local mesh size in structured environ¬ 
ments [iniiii]. In such cases the tracer diffusion becomes 
characterized by a non-uniform, position-dependent dif¬ 
fusivity D(x). Similarly, spatially varying transport char¬ 
acteristics are ubiquitous in contaminant dispersion in 
subsurface water aquifers [T2] . 

In the field of stochastic dynamics anomalous diffusion 
in spatially random media, disordered energy landscapes. 
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weakly chaotic systems, and dynamic maps received con- 
siderable attention [TSHSO] . More specifically, anomalous 
diffusion due to micro-domains was investigated [2T] and 
the influence of environmental Gaussian noise on diffus¬ 
ive particle trajectories in disordered systems was studied 
[22] . Moreover, deviations from normal diffusion due to 
quenched and annealed disorder of the medium diffusiv¬ 
ity received renewed interest [23J |21] . In such studies one 
is mainly interested in the quantitative behavior of the 
particle MSD ([^ as well as the ergodic properties of the 
system: is the information from time averages of physical 
observables typically garnered as time series by modern 
particle tracking assays equivalent to those of the corres¬ 
ponding ensemble averages known from the theoretical 
models? It turns out that a large variety of anomal¬ 
ous diffusion processes involve weak ergodicity breaking 
[3 UMin], the disparity between (long) time averages 
and ensemble averages of physical observables such as 
the MSD, and that in those cases the Khintchine theorem 
needs to be substituted by generalized versions [301 EH. 

Here, we study the dynamics and the ergodic prop¬ 
erties of heterogeneous diffusion processes (HDPs) with 
position dependent diffusivity D(x), in the presence of 
piece-wise deterministic quenched and annealed disorder. 
More specifically, we generalize HDPs with power-law dif¬ 
fusivity 

D{x)=Do\xr, ( 2 ) 

for which the anomalous diffusion exponent of the MSD 
assumes the form [53H5B] 


The physical dimension of the coefficient Dq in Eq. 0. is 
[Do] = cm^““°sec“^. The exponent 0 designates sub¬ 
diffusion for Qfo < 0 and superdiffusion for 0 < ao [32U36] . 
The profiles of the diffusivity for these cases are shown 
in Fig. 0 ,b. HDPs are weakly non-ergodic and ageing. 
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Figure 1: Particle diffusivity for heterogeneous diffusion pro¬ 
cesses: (a) and (b) for, respectively, ao = —2 and ao = 1- 
Panels (c) and (d) show the diffusivity of HDPs with quenched 
disorder for the same values of ao and for the parameters 
a'^ = 0.25 and 25a; = 2 (see text for details). Slight variations 
of the diffusivity in panels (c) and (d) indicate the external 
noise superimposed in the process. 


that is, their dynamics depends explicitly on the time 
gap between original initiation of the system and start of 
the measurement . We note that the ageing prop¬ 

erties of HDPs [5S] embodied in the ensemble and time 
averaged MSDs are in fact similar to those of subdiffusive 
continuous time random walks m and scaled Brownian 
motion [35]. 

In the following we unravel how the additional dis¬ 
order in the system modifies the diffusive and ergodic 
properties of HDPs. We compute the scaling laws for 
the ensemble and time averaged MSDs as well as the 
amplitude spread of individual realizations of the process. 
The article is structured as follows. In Sec. H we define 
the observables, that we will analyze. Sec. HI specifies 
the model and its implementation in the simulations. In 
Sec. jTVjwe then study HDPs with annealed disorder, fol¬ 
lowed by the scenario with quenched disorder in Sec. jV] 
Sec. VI concludes this work. 


II. OBSERVABLES 


The central quantity in the study of stochastic pro¬ 
cesses is the ensemble averaged MSD 


is the time averaged MSD 


52(A) 


1 

T-AJo 


x{t + A) 


r:{t) 


1 2 


dt. 


(5) 


Here A is the lag time and T is the measurement time 
(length) of the trajectory x{t) [3ll71|35j. Often, also the 
additional average 


1 ^ 

(^) = 


( 6 ) 


of the time averaged MSD over N individual trajectories 
is taken mEIlES]. A process is called ergodic when we 
observe the equality 


(x2(A))= lim 52(A). 
A/T-J-O 


(7) 


Examples for ergodic processes are Brownian motion 
as well as anomalous diffusion processes with 
MSD 0 given by random walks on fractals |40j and pro¬ 
cesses driven by fractional Gaussian noise [111133] . Once a 
process is non-stationary, the equality Q is violated, the 
phenomenon of weak ergodicity breaking [TlllSHinj- A 
whole range of anomalous diffusion processes with power- 
law MSD Q belongs to this class and specifically exhibits 
the linear lag time dependence 

(^) - (8) 

of the time averaged MSD [7]. As examples we mention 
continuous time random walk processes with scale free 
distributions of waiting times [3 dSHlSl EH, correlated 
continuous time random walks [U], as well as diffusion 
processes with space [321136] and time [321 |38l HU [46] 
dependent diffusion coefficients and their combinations 
m- We also mention ultraslow diffusion processes with 
a logarithmic form for {x‘^{t)) and linear lag time de¬ 
pendence Q of the time averaged MSD [48| as well as 
the ultraweak ergodicity breaking of superdiffusive Levy 
walks [13] . 

For finite measurement time even ergodic processes ex¬ 
hibit a statistical scatter of the amplitude of time aver¬ 
aged observables. This irreproducibility for the case of 
the time averaged MSD 52(A) can be quantified in terms 
of the distribution as function of the dimensionless 
variable [71 [33H35] 





(9) 


/ OO 

x'^P{x,t)dx, (4) 

-OO 

calculated as the spatial average of over the probab¬ 
ility density function P(x,t) to find the particle at po¬ 
sition X at given time t [39] . However, when individual 
time series x{t) of the particle position are measured in 
experiments or simulations, the typical quantity studied 


The variance of 0(^) is quantified in terms of the ergodi¬ 
city breaking parameter [71 l 26 ll 28 ] 

EB(A) = (^^(A)) - ^ (e) - 1. (10) 

For Brownian motion the behavior of the ergodicity 
breaking parameter at A/T —>■ 0 is 

4A 

EBbm(A) = —. (II) 
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Continuous time random walk processes with scale free 
waiting time distribution have a finite value for EB even 
in the limit A/T = 0 [5S], similar to HDPs [551 - I5S] . while 
for scaled Brownian motion the ergodicity breaking para¬ 
meter approaches zero in this limit [351135] . 

For reference in what follows we also mention that the 
probability density function of HDPs obeys has the ex¬ 
ponential form |33j 


I |-2/ao 

P{x,t) = exp 

^ ' ^/47r^ 


\2—ao 


{2-aoYD^t 


( 12 ) 


which is a stretched (compressed) Gaussian for super- 
diffusive (subdiffusive) HDPs with 0 < ao < 2 (ao < 0). 
Note that, respectively, the shape (12) has a distinct cusp 
at the origin or is bimodal with Pj^t) = 0 [55] . 


III. MODEL 

We employ the same tested stochastic algorithm for the 
Markovian HDPs as developed in Refs. |55H55] . based on 
the one-dimensional Langevin equation for the particle 
displacement x{t) with the position dependent diffusivity 
D{x), 


dx{t) 

dt 


= ^2D{x) X C(t). 


(13) 


The process is driven by the white Gaussian noise C(t) 
with covariance (C(^)C(^O) = ~ mean 

(C(t)) = 0. We interpret Eq. (13) in the Stratonovich 


sense leading to the following implicit mid-point iterative 
scheme: at step i -|- 1 the particle position is 




f ^i+1 “1“ 


12D 




x{yi+i-yi), (14) 


where the increments (y^+i — yi) of the Wiener process 
represent a (^-correlated Gaussian noise with unit vari¬ 
ance and zero mean. Unit time intervals separate con¬ 
secutive iteration steps. Below we simulate three values 
for the exponent ag, corresponding to j3 = 1/2 (sub- 
diffusive MSD), /? = 0 (Brownian motion), and /3 = 2 
(superdiffusive MSD). For standard HDPs these cases 
were analyzed by us in Refs. [33H36] . To avoid divergen¬ 
cies of the particle motion we regularize the diffusivity at 
x = 0 by addition of a small constant, namely D{x) = 
Do(|x|“ -I- Doff) where Doff = 10“^ and Dg = 10“^ for 
all results shown below. This choice does not affect the 
quality of the studied scaling laws [55] . 

To examine the effect of additional noise due to the 
environment we implement a Gaussian distribution of the 
scaling exponent of the diffusivity with the mean ag. 


p{a) 


1 


exp 


(a-ag)2\ 

2a2 J ■ 


(15) 


Generally, the distribution p{a) may be asymmetric, but 
we restrict our discussion to symmetric forms. We con¬ 
sider two versions of this additional disorder correspond¬ 
ing to the annealed and quenched limits for the variation 
of a. In the annealed case of noisy HDPs, the properties 
of the environment change rapidly in time compared to 
time scales of the particle motion. Physically, such noise 
may be due to the imprecision of the experimental setup 
or because of additional thermal agitation in the system. 
The diffusing particle thus visits regions in space with 
different local exponents a. In this scheme the particle 
diffusivity at position x fluctuates in time, and the value 
of the diffusivity will be different each time the particle 
revisits the same position x. In this annealed case large 
diffusivity variations occur in the entire space. 

In superdiffusive HDPs distant particle excursions take 
place due to the growth of D(x) away from the origin 
and the associated acceleration of the motion, while for 
subdiffusive HDPs the walker is increasingly trapped in 
the low-diffusivity regions at larger value of the position 
\x\ [33ll36] . With increasing strength cr^ of the annealed 
noise given by the distribution (151 the excursions of the 
particles in both superdiffusive and subdiffusive cases be¬ 
come more erratic as time evolves. The time interval 5t 
during which the walker has a given HDP exponent ai 
obviously affects the properties of noisy HDPs. These 
time spans St are here taken to be uniformly distributed. 
To simulate annealed noisy HDPs we use Eq. (15) with 
varying cr^. The particle performs jumps with a given 
scaling exponent for the time interval St, after which a 
new exponent is chosen from the distribution (15), and 
so on. The particle displacement Xi during the time span 
+ St) with HDP exponent Ui is the starting condi¬ 
tion for the next time interval. Shorter St intervals imply 
more erratic motion, as shown below. 

For noisy HDPs in the presence of quenched disorder, 
the profile of the particle diffusivity is hard-wired into the 
environment. We choose a static periodic arrangement of 
domains as shown in Fig. [^,d. In each domain the ex¬ 
ponent a is drawn from p{a) and the particle performs 
a regular HDP. The midpoint of each domain is chosen 
as the origin in the local HDP coordinate system, that 
is, locally the functional shape of D(x) is centered and 
decays or increases with the local scaling exponent a, as 
exemplified in Fig. [^^id. The period Sx for the strati¬ 
fication of the environment plays the role of a switch¬ 
ing mechanism affecting the system dynamics. At the 
boundary of the domains the diffusivity and its derivat¬ 
ive in general acquire jumps. Physically, the latter occurs 
in the presence of some walls, cages, etc. 

We simulate quenched noisy HDPs as follows. The en¬ 
tire space is stratified into domains of width 2Sx, and 
the local HDP exponent is chosen from the distribution 
(15). The length Sx is a vital parameter of quenched 
noisy HDPs. The particle performs an HDP random 
walk in each space domain with D(x, a) and it hops 
to a neighboring domain once the domain boundary is 
reached. The centers of the domains are computed from 
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the particle position Xi as 


Xf.,i = 2{dx)mt 


25x 


-sign[a;i](fe), 


(16) 


see Fig. [^,d. Here int[a:] denotes the integer part of the 
argument, and an additional Sx shift is used for conveni¬ 
ence. The starting position of the particle is near the 
center of the first domain, at a;(0) = 0.1 -|- Xc,i- The sub¬ 
sequent position Xi+i is evaluated from Xi with the local 
exponent Ui according to Eq. (14), that is, 


— \ 2,^0 


Xi-\-\ Xi 


Doff 


x(2/*+i - y^)■ 


(17) 


We vary the width of p{a) and the mean value of the 
scaling exponent a^. Shorter periodicities 5x are equi¬ 
valent to stronger external noise, as shown below. We 
note here that for subdiffusive HDPs the centers of the 
domains Xc^i correspond to the regions of maximal diffus- 
ivity, while for superdiffusive HDPs these are the spots 
of the lowest diffusivity [S3 IMl • 


IV. NOISY HDPS WITH ANNEALED 
DISORDER 


For even larger noise strength cr^, the behavior of the 
time averaged MSD and the ergodic properties are dom¬ 
inated by extreme events, that is, by single or few traject¬ 
ories in the data set with the largest exponent (s) yielding 
extremely distant particle excursions. With an increasing 
width of the a distribution p(a), the spread of the time 
averaged MSD grows, as well, as evidenced in Fig. 
Similarly, for such large values of the noise strength cr^ 
the value of the ergodicity breaking parameter becomes 
proportional to the number N of recorded traces, wit¬ 
nessing the dominance of single traces, each having the 
potential to be more extreme than the others, compare 
Fig. [3^. We refer the reader to Ref. [3S], in which the 
critical properties of HDPs and the effects of the number 
of traces are analyzed in the limit ag — > 2. 

For narrow distributions pia) the spread <(>(^) of indi¬ 
vidual 6“^ traces is symmetric at short lag times A, de¬ 
veloping a tail at larger lag times A. This behavior can 
be rationalized in terms of a generalized Gamma distri¬ 
bution (see Ref. [33 )■ The general features of are 
shown in Fig. in terms of the higher moments of this 
distribution. These are the skewness 


* ^ ~ N ‘ \ 3/2 

and the kurtosis 


(18) 


A. Noisy Brownian motion, ao = 0 


For ag = 0 and a small value cr^ of the additional noise, 
as expected, we observe small discrepancies from the ca¬ 
nonical Brownian motion, as evidenced in Fig. [2)o. The 
behavior is ergodic, and the ergodicity breaking para¬ 
meter follows the known behavior (11) for Brownian mo¬ 
tion, see Fig. [^ Most importantly, the ensemble aver¬ 
aged MSD equals the time averaged MSD, apart from 
very short lag times at which the relaxation from the ini¬ 
tial value a;(0) = Xq occurs (compare Ref. [35] for more 
details). At longer lag times, the deteriorating statistics 
of the traces give rise to the typical cone-like scatter. 

As the noise strength cr^, the variance of the a distri¬ 
bution p(a) is increased, see, for instance, in the panel 
for the noise strength a — 0.5 in Fig. a more pro¬ 
nounced scatter of the 6“^ traces emerges and, import¬ 
antly, persists in the limit A/T —> 0. The occurrence 
of progressively more distant particle excursions caused 
by superdiffusive traces with a > 0 gives rise to a larger 
spread of the amplitude scatter quantified by the distri¬ 
bution The value of grows somewhat faster 


than the ensemble MSD ([^ due to these outliers, giving 
rise to larger values of the ergodicity breaking parameter 
EB (not shown). The time averaged MSD S'^{A) scales 
linearly with the lag time A, and, as they should, in the 
limit A —>■ T the time averaged MSD settles back to the 
ensemble averaged MSD, due to the pole in the definition 
([^ of the time average. 


K{0 = 




(19) 


which complement the variance of described by the 
ergodicity breaking parameter (10). In Eig.[^ we also ob¬ 
serve that for small noise strengths tr^ the value EB(A = 
1) for noisy Brownian motion approaches EBbm(A = 1) 
given by Eq. as expected. The values of the er¬ 

godicity breaking parameter grow with A, indicative of a 
bigger spread of the value 6“^ of individual traces (green 
points in Fig. [^). 


B. Subdiffusive noisy HDP, oq = —2 

For the subdiffusive case the time evolution of the en¬ 
semble and time averaged MSDs is illustrated in Fig. 
for different noise strengths of the a distribution. We 
observe that for the subdiffusive value ag = —2 the same 
magnitude of the a variation causes a much weaker effect 
as compared to the Brownian (ag = 0) or superdiffusive 
(ag = 1) situations. The scatter of remains nearly 
insensitive to the lag time A, similar to canonical HDPs 
[531135] . The scaling of the ensemble averaged MSD also 
agrees with that for HDPs [33] ■ It is reached after less 
than a dozen of steps during which the relaxation of the 
initial condition occurs, compare Refs. [53JI35]- The scal¬ 
ing of the time averaged MSD (S'^) remains linear and 
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Figure 2: Ensemble and time averaged MSDs and ^(5^(A)y (thick blue curves) as well as individual time averaged MSDs 

(red curves) for annealed noisy HDPs. Parameters: the values of ao and its variance are indicated in the plots, the trace 
length is T = 10'*, and the number of sampled traces is N = 10^. The initial position is xo = x{t = 0) = 0.1. The top panels 
correspond to the noisy subdiffusive case, the middle panel represents noisy Brownian motion, and the bottom panels are the 
case of superdiffusive noisy HDPs. The asymptotes 0 and ® for the ensemble and time averaged MSDs of standard HDPs 
are shown as the dashed curves. The Brownian asymptote (x^t)) = 2Dot is the dashed-dotted line. 
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Figure 3: Second, third, and fourth order moments of the amplitude scatter distribution (^(^) for annealed noisy HDPs, 
computed for the parameters of Fig. Large, medium, and small symbols correspond to lag times A = 1, 10, and 100, 
respectively. The dotted line in panel (a) at large is the ergodicity breaking parameter EB ^ N = 10^, indicative of the 
single trace dominance in this case, see text. The dotted lines in panel (a) for small noise strength stands for EBbm(A) 
given by Eq. (111. 


nearly unaffected by changes of a^. The long time scal¬ 
ing of the MSD is also weakly sensitive to in the range 
considered here. 

Physically, for the subdiffusive case the spread of 
should be > Oq to give rise to fast particle excursions 
(outliers). Thus, much larger values are required to 
disturb the spread of for strongly subdiffusive noisy 
HDPs as compared to superdiffusive noisy HDPs shown 
in Fig. [^. This is our first important conclusion. 

We rationalize the effects of the a spread further in 
terms of the width and higher moments of the amplitude 
scatter distribution The results for ao = 0, sub- 

and superdiffusive annealed noisy HDPs are shown in 
Fig.[3]a. We observe that all moments are typically smal¬ 
ler for the subdiffusive case reflecting a less pronounced 
and asymmetric spread of . The skewness of Brownian 


motion —?► 0) tends to vanish, as it should, while 

for sub- and superdiffusive noisy HDPs it attains finite 
values at cr^ —)■ 0 (Fig. ^). This is due to the inher¬ 
ent asymmetry of the cj)(^) scatter even at cr^ —>■ 0: it 
features a tail at large ^ values, a maximum at inter¬ 
mediate and vanishes at ^ 0 [33]. Both skewness 

S'(^) and kurtosis K{^) grow dramatically with cr^ for all 
values of ao, as demonstrated in Figs. [3]d,c. We checked 
that for (T^ —0 the value of the ergodicity breaking para¬ 
meter in the limit A/T 1 approaches that for standard 
HDPs |33| , as expected, while for a broad distribution of 
a values the ergodicity breaking parameter increases and 
eventually approaches the number of traces TV in the data 
set (single-trace domination), Fig. The value of the 
ergodicity breaking parameter EB for oq < 0 is nearly 
unaffected by variations of tr^ over a wide range, see the 







7 


red symbols in Fig. This reflects the minor change 
in the spread of single traces 5'^ when cr^ is varied, see 

Fig. 1^. 


C. SuperdifFusive noisy HDPs, ao = 1 

The ensemble and time averaged MSDs of superdiffus- 
ive noisy HDPs with ao = 1 are shown in Fig. For 
small noise strengths their scaling agrees with the res¬ 
ults for standard HDPs, Eqs. Q and Q. With increas¬ 
ing noise strength the time averaged MSD traces (5^ 
grow dramatically, and for moderate and large lag times 
A the time averaged MSD deviates progressively from 
the HDP scaling, that is ballistic for uq = 1 (Fi g. [2] :). 
The scatter of the individual time averaged MSDs 5“^ be¬ 
comes progressively larger and asymmetric as the width 
of p{a) increases. The amplitude of the time averaged 
MSD traces for large values of grows significantly 
above the asymptote for undisturbed HDPs due to single 
trajectory domination. Therefore, the moments of the 
scatter distribution increase, see the blue symbols 
in Fig. For the later parts of the trajectories, the en¬ 
semble averaged MSD increases very fast (see the right 
panel in Fig. [^) to meet the value of in the limit 
A = T. For superdiffusive HDPs the moments of 
are larger than those for subdiffusive noisy HDPs with 
the same tr^, compare the red and blue symbols in Fig.[^ 


V. NOISY HDPS WITH QUENCHED 
DISORDER 


We now turn to the situation of quenched disorder in 
a stratified environment, in which evenly sized domains 
of width 6x have a diffusivity of the form Q, centered 
within the domain, whose a value is noisy and with dis¬ 
tribution (15). In this quenched scenario, the particle 
experiences the same value of a each time it revisits a 
given domain. The situation is illustrated in Fig. [^,d. 


A. Noisy Brownian motion, oq = 0 


For quenched, noisy Brownian motion we observe that 
for small noise strength the behavior, as expected, 
is very close to standard Brownian motion (not shown). 
For a large value of cr^, the spread of the amplitude of 
individual time averaged MSDs (5^ (A) is non-negligible 
even at short lag times A, as shown in Fig. H- This 
spread is more pronounced for larger periodicities 5x of 
the stratified medium. For small cr^ the ensemble av¬ 
eraged MSD ((x — Xc,i)^) computed with respect to the 
center of the starting domain and the time averaged MSD 
(thick blue lines) almost coincide for all lag times A (not 
shown). Concurrently, the ergodicity breaking parameter 
follows the Brownian asymptote (11), as shown by the 


green symbols in Fig. For larger values of the noise 
strength the ergodicity breaking parameter deviates 
pronouncedly from Eq. ( |11[ ) at short lag times A, indic¬ 
ating the occurrence of weak ergodicity breaking, along 
with the disparity ^ (^^)> witnessed by Eig.j^. 
This inequality is particularly pronounced for larger val¬ 
ues of the noise strength and large periodicity fe, see 
the changes for varying bx in Eig. Eor wider a. dis¬ 
tributions via) the ensemble averaged MSD starts close 
to that of the asymptote for standard Brownian motion, 
while at later times there occurs a crossover to the curve 
for the time averaged MSD (left panel, Eig. |^.) This be¬ 
havior is also typical for sub- and superdiffusive quenched 
noisy HDPs, see below. For = 1 this transition occurs 
after ~ 10^ time steps and becomes less pronounced for 
smaller periodicities bx of the medium (Fig. §)• 


B. Superdiffusive noisy HDPs, oq = 1 


In standard superdiffusive HDPs there exists a finite 
probability of particle trapping in regions of low diffus- 
ivity near the origin, as witnessed by the cusp around 
X = 0 of the probability density function (12) [33]. For 
noisy HDPs we find that for large values of the domain 
size bx and small noise strengths the particle prefer¬ 
entially stays in the domain, in which it was seeded, and 
the resulting ensemble averaged MSD is close to that of 
the standard HDPs [531 Here we again computed the 
MSD with respect to the center Xc,i of the seed domain in 
the form x{t = 0) — Xc,i = 0.1. The time averaged MSD 
is equally close to the asymptote (|^ of the normal HDP. 
Ensemble and time averaged MSDs converge at long lag 
times A —T, note that the ensemble averaged MSD 
here is below the time averaged MSD, as evidenced by 
Eigs. |§i,b. 


We start with a narrow spread of a in the spatial do¬ 
mains corresponding to u = 0.03. In this case we find 
that with decreasing domain size bx the amplitude scat¬ 
ter of individual time averaged MSDs shrinks and the 


amplitude of the trajectory mean drops sub¬ 

stantially (Fig. 1^). The reason is that for a small do¬ 
main size there are almost no regions of fast diffusivity. 
For small values of bx the ensemble and time averaged 
MSDs converge and drop below the Brownian asymp¬ 
tote, see the dashed-dotted line in the right graph in 
Fig. [^. In such cases of smaller domain size the ergodi¬ 
city breaking parameter attains relatively small values, 
as shown in Fig. [^, indicating a more ergodic behavior. 
This effect of the noise is similar to that for noisy CTRWs 
[22] . As bx increases, the ergodicity breaking parameter 
approaches values close to those of the standard HDP, 
EB(A = 1) « 0.34 for ao = -2 and EB(A = I) « 1.1 for 
ao = 1, with T = 10'^ [33l[35|. This is indicated by the 
dashed-dotted lines in Fig. Hi- Thus, frequent hopping 
events between individual domains destroys the charac¬ 
teristic of the noise-free HDP scaling and causes the dif- 











Figure 4: Ensemble and time averaged MSDs and amplitude scatter of individual traces <5^ for noisy HDPs with quenched 
disorder. The values of a, a, and domain size Sx are indicated in the plots. The panels (a) and (b) are for subdiffusive noisy 
HDPs, panel (c) stands for noisy Brownian motion, and panels (d) and (e) represent superdiffusive noisy HDPs. The MSD is 
computed with respect to the position of the center of the first domain, {{x{t) — Xc,i)'^)- Parameters: T = 10"^, N = 10^, and 
5x values are the same in each column. The notations for the curves and asymptotes are the same as in Fig. 
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Figure 5: Ergodicity breaking parameter of noisy HDPs with quenched disorder. The parameters are the same as in Fig. 
the values for a and cr being indicated in the plots. The black dotted lines represent the ergodicity breaking parameter for the 
standard HDPs. The large, medium, and small symbols correspond to the lag times A = 1, 10, and 100, respectively. 


fusion to be more ergodic. This is our second important 
conclusion. 

For larger a values the MSD stops following the HDP 
scaling law Q and instead two nearly Brownian regimes 
are detected for short and long diffusion times, see the 
left panel in Fig. |^. Similar to noisy CTRWs [71 
for noisy HDPs we observe a superposition of anomal¬ 
ous scaling for the MSD inherent to HDPs with the lin¬ 
ear MSD increase due to particle jumping between the 
stratified domains. The latter term contributes stronger 
for smaller 5x values: after a given number of steps T 
performed the particle visits more D{x) domains and its 
diffusion on the length scale Sx becomes effectively 
more normal and ergodic. 

The time averaged MSD is an approximately linear 
function of the lag time A. For smaller domain size Sx we 
observe a more confined amplitude spread of the time av¬ 
eraged MSD traces S^, see the evolution from left to right 
in Figs. |^,b. A similar behavior occurs for subdiffus- 
ive noisy HDPs, as demonstrated in Figs. |^,e consistent 
with smaller values of the ergodicity breaking parameter. 
This is our third main result. For superdiffusive noisy 
HDPs, given large domain sizes Sx, we observe more dis¬ 
tant particle excursions and thu s a br oader amplitude 
spread of individual time traces <5^(A), particularly for 
large values of the noise strength (Figs.|^): at larger 
we correspondingly obtain larger values of the ergodi¬ 
city breaking parameter, compare panels (a) and (b) in 
Fig. [5] 


We find that the distribution i^(^) of the amplitude 
scatter features a skewed form, which is characterized by 
its second, third, and fourth moments corresponding to 
the ergodicity breaking parameter, the skewness S, and 
the kurtosis K, respectively. For larger values of S(^) 
and Ar(^) grow with the domain size Sx and are more ir¬ 
regular than the distribution itself, due to worsening 
statistics for higher order moments (not shown). Note 
that for short lag times A/T <C 1 the ergodicity break¬ 
ing parameter for large domain sizes Sx approaches the 
values of the corresponding normal HDPs [33], compare 
Fig. [§i. At small domain size Sx the non-ergodic prop¬ 
erties of the standard HDPs are masked by the noise in 
the stratified spatial domains. 


For superdiffusive HDPs the particles tend to localize 
in the center of each domain, while for subdiffusive values 
ao < 0 they tend to spread towards regions of low dif- 
fusivity near the domain borders. In the long time limit 
the particles spread over many domains, establishing the 
shape of the probability density function P{x,t) presen¬ 
ted in Fig. 1^ The local minima and maxima of P{x,t) 
correspond to the regions of low and fast diffusivity D{x), 
respectively, see Figs.[^,d and|^ For relatively large do¬ 
main size Sx the probability density function of the noisy 
HDPs becomes dominated by the contribution from the 
seed domain. The spreading of particles over superdif¬ 
fusive HDP domains in the long time limit is symmetric 
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Figure 6: Probability density functions of noisy HDPs with quenched disorder for noise strength = l and varying domain size 
5x. The other parameters are the same as in Fig. Panel (a) represents superdiffusive noisy HDPs and panel (b) stands for 
subdiffusive noisy HDPs. The shift of the pe ak p ositions with respect to those of the standard HDPs [33] is due to the shift 
(dx) of the domain center positions, see Eq. ( |16| ). The dashed curves represent the probability density functions of standard 
HDPs, given by Eq. (12 1 , while the dashed-dotted curves are result (20l. 


and nearly Gaussian, 


P{x,t) 


1 

, - exp 

y/47r£)efft 


V 4Deff^ / 


( 20 ) 


with the effective diffusivity Heff- The mean particle dis¬ 
placement with respect to the center of the seed domain 
vanishes, {x(t —?> oo)) —?► 0. To compute Heff analytically 
a homogenization procedure and generic concepts of dif¬ 
fusion in random and highly heterogeneous media would 
need to be applied [50] . 


C. Subdiffusive noisy HDPs, oq = —2 

Subdiffusive noisy HDPs in the quenched scenario 
share a number of trends with the above descriptions of 


the cases a = 0 and a = 1. In particular, as the domain 
size 5x decreases, the cimplitude spread of individual time 
averaged MSD traces <5^ decreases (Fig.|^,e). Because of 
the sublinear scaling of the ensemble MSD of the normal 
subdiffusive HDPs (oq < 0) the ensemble averaged MSD 

approaches the time averaged MSD ^(5^(A)^ from above. 
Moreover, the scaling of the ensemble averaged MSD of 
subdiffusive noisy HDPs with quenched disorder turns 
from subdiffusive to Brownian as the domain size 5x de¬ 
creases. The physical reason for this crossover behavior 
is the random character of hops between domains with 
a varying local exponent a. We find that, similarly to 
superdiffusive noisy HDPs, the ensemble averaged MSD 
initially follows the scaling 0 of normal HDPs while at 
later times a nearly linear scaling is observed. For smal¬ 
ler periodicities 5x the linear scaling becomes dominant, 
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as demonstrated in Fig. from left to right. 

The probability density function of quenched noisy 
HDPs in the long time limit is a combination of the super¬ 
imposed local probability densities of the standard HDP. 
For large periodicities 6x the probability density function 
is again dominated by the contribution from the seed do¬ 
main, as can be seen in the right panel of Fig. Similar 
to superdiffusive noisy HDPs we find that the time aver¬ 
aged MSD is linear in the lag time, ^ A, while 

the amplitude spread of individual time averaged MSDs 
grows with the noise strength and becomes diminished 
for smaller medium periodicities 5x. We also see that for 
subdiffusive noisy HDPs the saturation of the ergodicity 
breaking parameter to the values of normal HDPs occurs 
at much smaller values of 6x as compared to superdiffus- 
ive noisy HDPs (Fig. [^). 

VI. CONCLUSIONS 

We studied a stochastic process based on a combina¬ 
tion of heterogeneous diffusion processes with multiplic¬ 
ative noise and additional disorder of the environment, 
distinguishing annealed and quenched scenarios. The en¬ 
vironment was assumed to be structured into periodic 
domains of given periodicity. We investigated the diffus¬ 
ive and ergodic properties of these noisy heterogeneous 
diffusion processes. The superposition of the additional 
stochasticity onto the standard HDP with its determ¬ 
inistic variation of the diffusivity revealed a variety of 
new features, the scaling relations for the ensemble and 
time averaged MSD of the noisy HDPs being dramatic¬ 
ally altered as compared to the normal HDP behavior. 

For annealed disorder, the scaling exponent a of the 
diffusivity profile switches in time and the gradient field 
of the particle diffusivity has a single origin at a; = 0. 
We demonstrated how the Gaussian spread p{a) of the 
scaling exponent gives rise to a strongly asymmetric scat¬ 
ter of individual time averaged MSD traces. Rapidly 
switching diffusivity profiles in such an annealed environ¬ 
ment cause transient particle trapping in low-diffusivity 
regions. For superdiffusive motion the effects of the a 
spread are more pronounced. In the case of a quenched 
environment, a spatially stratified medium is modeled in 
terms of domains of width 25x with a normal distribution 
of the local HDP exponent. Upon particle diffusion, the 
averaging is thus performed over ensembles of particle 
trajectories generated for different spatial distributions 
of the scaling exponents a in the domains. One of the 
key findings is that for small periodicity 5x the sub- and 
superdiffusive scaling of normal HDPs cross over to a lin¬ 
ear growth of the ensemble averaged MSD as function of 
time. External noise thus progressively masks the stat¬ 
istics of the underlying HDP. 

What could be the physical phenomena captured by 


the noisy HDP discussed here? From a biological per¬ 
spective, the diffusion of small molecules in assemblies 
of non-identical, interconnected cells is a relevant ex¬ 
ample. The cell-to-cell variations of the diffusivity are 
inherent to biological tissues, while every individual cell 
features a space dependent diffusivity in its cytoplasm 
[5]. At cell-to-cell boundaries the diffusivity likely var¬ 
ies with a jump, as captured by our stratified model of 
the quenched disorder, with possibly discontinuous dif¬ 
fusivity across the system. We note that heterogeneous 
diffusivities can, for instance, play a role in the formation 
of gradients of morphogen molecules in a developing cell 
tissue m, a process known to involve features of anom¬ 
alous diffusion. It also features a division of fluxes of the 
molecules into fluxes through cells, across the outer cell 
membranes, and transport in extracellular spaces [52) . 
Heterogeneous diffusion of water molecules in brain tis¬ 
sues [53] and strongly heterogeneous structures of cardiac 
muscle tissue with nontrivial cell-cell coupling [^ could 
be another example. Similarly, the domains in the noisy 
HDP could represent internal compartments in a single 
cell. The quenched case would correspond to static en¬ 
vironments whereas the annealed scenario would stand 
for environments, which change rapidly compared to the 
typical crossing times between domains. 

Our results for noisy HDPs could also be useful for the 
description of nano-objects trapped in dynamical tem¬ 
perature fields [5S] and of particles in strong temper¬ 
ature gradients |S^. Another field of relevance is the 
tracer diffusion in heterogeneous assemblies of distrib¬ 
uted obstacles m mimicking features of the cell cyto¬ 
plasm |S] and diffusion on chemically and mesoscopically 
periodically patterned solid-liquid interfaces jSS]. On a 
macroscopic scale, water diffusion in subsurface hydro¬ 
logy applications is to be mentioned as well as tracer 
motion in porous heterogeneous media [59] . For the latter 
there likely exists a distance-dependent diffusivity within 
each pore constructing a network governing the diffusion 
of water and contaminants in soil specimen [12] . Finally, 
in statistical models of financial stock price variations (50] 
the terms stochastic versus correlated volatility widely 
occur, representing the diffusivity in random walk mod¬ 
els m- Some patterns of correlated or clustered volat¬ 
ility observed in financial data thus correspond to a sys¬ 
tematically varying diffusivity in our model of quenched 
noisy HDPs. Some repeats of non-Brownian up-and- 
down trends in stock price fluctuations m can thus be 
considered as HDPs repeatedly occurring in time. 
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